/***************************************************************************//**
 * @file FDM/FDM_Atmosphere.cpp
 * @author  Marek M. Cel <marekcel@mscsim.org>
 *
 * @section LICENSE
 *
 * Copyright (C) 2013 Marek M. Cel
 *
 * This file is part of MScSim. You can redistribute and modify it under
 * the terms of GNU General Public License as published by the Free Software
 * Foundation; either version 3 of the License, or (at your option) any later
 * version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
 * for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.
 * 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 * Further information about the GNU General Public License can also be found
 * on the world wide web at http://www.gnu.org.
 ******************************************************************************/
#ifndef FDM_ATMOSPHERE_CPP
#define FDM_ATMOSPHERE_CPP
#endif

////////////////////////////////////////////////////////////////////////////////

#include <cmath>
#include <iostream>
#include <stdlib.h>

#include <FDM/FDM_Atmosphere.h>

////////////////////////////////////////////////////////////////////////////////

using namespace std;

////////////////////////////////////////////////////////////////////////////////

const double FDM_Atmosphere::m_G = 9.80665;
const double FDM_Atmosphere::m_R = 8.31432e3;
const double FDM_Atmosphere::m_S = 110.4;
const double FDM_Atmosphere::m_B = 1.458e-6;

////////////////////////////////////////////////////////////////////////////////

FDM_Atmosphere::FDM_Atmosphere()
{
    // standard sea level conditions
    m_tempK_0 = 288.15f;
    m_press_0 = 101325.0f;
    
    m_tempRate = 6.5e-3f;
    
    m_tempK = 288.15f;
    m_tempC = 15.0f;
    m_press = 101325.0f;
    m_dens  = 1.225f;
    m_sos   = 340.29f;
    m_visc  = 1.79e-5f;
}

////////////////////////////////////////////////////////////////////////////////

FDM_Atmosphere::~FDM_Atmosphere() {}

////////////////////////////////////////////////////////////////////////////////

void FDM_Atmosphere::init( double tempK, double press, double alt )
{
    if ( alt < 11000.0 )
    {
        // temperature drop rate up to 11 000 m ASL
        // calculated due to variable surface temperature
        m_tempRate = ( tempK - 216.65 ) / ( 11000.0 - alt );

        m_tempK_0 = 11000.0 * m_tempRate + 216.65;
        m_press_0 = press / pow( 1.0 - m_tempRate * alt / m_tempK_0, 5.2561 );
        
        update( alt );
    }
    else
    {
        cerr << "Error in Atmosphere::init( float, float, float ): Atmosphere initial altitude over limit." << endl;
        exit( EXIT_FAILURE );
    }
}

////////////////////////////////////////////////////////////////////////////////

void FDM_Atmosphere::update( double alt )
{
    if ( alt <= 82000.0 )
    {
        if ( alt <= 11000.0 )
        {
            // temperature and pressure
            m_tempK = m_tempK_0 - m_tempRate * alt;
            m_press = m_press_0 * pow( 1.0 - m_tempRate * alt / m_tempK_0, 5.2561 );
        }
        else
        {
            // temperature
            if ( alt > 11000.0 && alt <= 20000.0 ) m_tempK = 216.65;                                // const
            if ( alt > 20000.0 && alt <= 32000.0 ) m_tempK = 216.65 + ( alt - 20000.0 ) * 1.0e-3;   // const rate +1.0 K/km
            if ( alt > 32000.0 && alt <= 47000.0 ) m_tempK = 228.65 + ( alt - 32000.0 ) * 2.8e-3;   // const rate +2.8 K/km
            if ( alt > 47000.0 && alt <= 51000.0 ) m_tempK = 270.65;                                // const
            if ( alt > 51000.0 && alt <= 71000.0 ) m_tempK = 270.65 - ( alt - 51000.0 ) * 2.8e-3;   // const rate -2.8 K/km
            if ( alt > 71000.0 && alt <= 82000.0 ) m_tempK = 214.65 - ( alt - 71000.0 ) * 2.0e-3;   // const rate -2.0 K/km
            
            // pressure
            m_press = m_press_0
                   * exp( -m_G * ( alt - 11000.0 ) / ( 287.05 * m_tempK ) )
                   * pow( 1.0 - 0.0065 * 11000.0 / m_tempK_0, 5.2561 );
        }
        
        m_tempC = m_tempK - 273.15;                                     // [deg C] temperature
        m_dens  = m_press / 287.05 / m_tempK;                           // [kg/m^3] density
        m_sos   = 331.3 * sqrt( 1.0 + ( m_tempK - 273.15 ) / 273.15 );  // [m/s] speed of sound
        m_visc  = m_B * pow( m_tempK, 3.0 / 2.0 )/( m_tempK + m_S );    // [Pa*s] dynamic viscosity
    }
    else
    {
        cerr << "Error in Atmosphere::update( float ): Atmosphere current altitude over limit." << endl;
        exit( EXIT_FAILURE );
    }
}
